Universality in sandpiles 
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We perform extensive numerical simulations of different versions of tfie sandpile model. We find 
that previous claims about universality classes are unfounded, since the method previously employed 
to analyze the data suffered a systematic bias. We identify the correct scaling behavior and conclude 
that sandpiles with stochastic and deterministic toppling rules belong to the same universality class. 
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Sandpile automata 0| are among the simplest mod- 
els to describe avalanche propagation, a phenomenon of 
upsurging experimental interest in a wide range of fields 
1^. In the stationary state, after suitable tuning of the 
driving fields Q, these models display critical behavior 
in the avalanche statistics. As for ordinary critical phe- 
nomena, it is possible to define a set of scaling exponents 
that characterize that large scale behavior of the system 

i- 

The precise identification of universality classes in 
sandpile models |^] is an unresolved issue. From a theo- 
retical standpoint, it would be unusual that small modi- 
fications in the dynamical rules of the model could lead 
to different universality classes. Real-space renormaliza- 
tion group calculations ||] suggest that different sandpile 
models, such as the Bak, Tang and Wiesenfeld (BTW) |^] 
and the Manna Q models, all belong to the same univer- 
sality class. This result is also confirmed by a recently- 
proposed field theory approach |6| , which shows that all 
sandpile models are described by the same effective 
field theory at the coarse grained level. Universality is 
also found between BTW (discrete) and Zhang Q (con- 
tinuous) models in the dynamical renormalization group 
calculations of Ref. ||^ . 

The results obtained by numerical simulations are 
unclear. Early large scale numerical simulations of 
the Manna @] and BTW models show that the 

avalanche distributions are described by the same ex- 
ponents for the power law decay and the scaling of the 
cutoffs. These results were questioned by Ben Hur and 
Biham who analyzed the scaling of conditional ex- 
pectation values of various quantities. They found 
significant differences in the exponents for the two models 
and therefore proposed a new classification of universality 
in sandpile models in which models with stochastic up- 
date rules, such as the Manna model, fall in a universality 
class different from that of Abelian models, such as the 
BTW fill . The method was later applied to the Zhang 
model that was declared "non- universal" This re- 

sults pose a puzzling problem, since they contradict all 
the existing theories and do not agree with the scaling 
predicted analyzing avalanche distributions . 



Here we present large scale numerical simulations of 
the BTW and Manna sandpile models, with the goal of 
settling the issue of universality. First we show that the 
method of conditional expectation values, introduced in 
Ref. Q and used in Ref. is systematically biased 
by non-universal corrections and does not provide indi- 
cations on universality classes. By removing the bias, we 
provide evidence that the BTW and Manna models are 
universal. We confirm this conclusion by data collapse 
and moment analysis of the distributions [p^ . 

Sandpile models are defined on d— dimensional hyper- 
cubic lattice. On each site i of the lattice we define an 
integer variable Zi which we call "energy" . At each time 
step an energy grain is added on a randomly chosen site 
{zi Zi + 1). When one of the sites reaches or exceeds 
a threshold Zc a "toppling" occurs: Zi = Zi — Zc and 
Zj ^ Zj + 1, where j represents the nearest neighbor sites 
of site i. In the BTW model Zc = 2d and each nearest 
neighbor receives a grain after the toppling of the site i. 
In the Manna model Zc = 2 and therefore only two ran- 
domly chosen neighboring sites receive a grain. A top- 
pling can induce nearest-neighbor sites to topple on their 
turn and so on, until all the lattice sites are below the 
critical threshold. This process is called an avalanche. A 
slow driving is usually imposed, so that grains are added 
only when all the sites are below the threshold. The 
model is conservative and energy is dissipated only at 
boundary sites Here, we perform numerical simula- 
tions of two-dimensional Manna and BTW models with 
open boundary conditions and conservative dynamics. 
The lattice size ranges from L — 128 to L = 2048 in 
both models. In each case, statistical distributions are 
obtained averaging over 10^ nonzero avalanches. 

Avalanches in sandpile models are usually character- 
ized by three variables: the number of topplings s, the 
area a affected by the avalanche, and the avalanche du- 
ration T. The probability distribution of each of these 
variables is usually described as a power law with a cutoff 



P{x) = X-^'^Gix/Xc), 



(1) 



where x = s,a, T. When the system size L goes to infin- 
ity the cutoff Xc diverges as Xc ^ L^^ ■ Under the finite 
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size scaling (FSS) assumption of Eq. (|l|), the set of expo- 
nents {tx, (3x} defines the universality class of the model. 

In two dimensions, an accurate numerical determina- 
tion of the power law exponents in Eq.(^ proved to be 
a difficult task [^| jl0|jl^Jl7|] , due to the large deviations 
at the lower and upper cutoffs. For this reason, Chris- 
tensen et al. in order to distinguish among univer- 
sality classes proposed a more refined numerical analysis 
based on the evaluation of the expectation value E{x\y) 
of the variable x restricted to all the avalanches with vari- 
able Y = y, where {X, Y} = {s, a, T} It is assumed 
that E{x\y) ^ y^'^y and the exponents ■jxy are used to dis- 
tinguish among universality classes |0| . These exponents 
satisfy the scaling relations -fxy = ly^ and -f^z = Ixylyz- 

If the conditional probability distribution p(a;|y) is suf- 
ficiently peaked, then ^^y is well-defined and to each 
value of the variable x we can unambiguously associate 



E{x\y)^ / dx-f{{x-y)/x*). 



a value of the variable y (i.e. 



y^^y). In partic- 



ular, the cutoff of the distributions should be related 
by the same exponents (i.e. Xc ^ yc^"), which implies 
Ixy = Px/fiy For instance, we have = f3s/2, since 
in two dimensions avalanches are compact for both the 
BTW ^ and Manna model Q, so that Pa = 2. The 
data collapse analysis shows the BTW and Manna model 
both share the same exponent f3s — 2.7 |Q, lOjl^ which 
implies 7^0 — 1.35. On the contrary, Refs. |ll]l4l] found 
7sa ^ 1-06 for the BTW model and -fsa = 1-24 for the 
Manna model, which would yield two different universal- 
ity classes for the two models. Less marked differences 
were also observed for the other exponents 7a;y |ll| , |l4| . 

In order to resolve this paradox, we return to the hy- 
pothesis underlying the use of conditional expectation 
values: p{x\y) must be symmetric and strongly peaked 
around the average value. We checked numerically that 
this assumption is not fulfilled: in the BTW model the 
distribution p{s\a) is maximum for s = a and decreases 
for s > a, with a characteristic value s* scaling as a'^^/^ 
(see Fig. 1). The distribution is not symmetric (see 
also Ref. U^), consistent with the constraint s > a 
(the avalanche area can not be greater than the num- 
ber of topplings). Similar considerations apply as well to 
other quantities (i.e. a > T, s > T) whose conditional 
probability distributions show asymmetry although less 
marked. 

To understand the effect of non-symmetric distribu- 
tions on conditional expectation values, consider a dis- 
tribution of the form 



p{x\y) = d{x - y)f{{x - y)/x*)/x* 



(2) 



where the characteristic value scales as x*{y) ^ y^"^" and 
9{x) is the step function . The factor 1/x* ensures nor- 
malization for any y 



dxp(x\y) = 1, (3) 



so that the conditional expectation value is given by 



(4) 



Performing the substitution z = x — y, we obtain 

z 

E{x\y)=y+ dz—p{z/x*)=y + Cy''-y, (5) 
Jo ^ 

where C is a non- universal constant. 

In the BTW model p{s\a) has the form of Eq. (||) as 
is shown by performing data collapse analysis (see inset 
of Fig. 1). Thus, we can easily subtract the hnear bias 
from the expectation value in order to obtain the cor- 
rect scaling behavior to be compared with that of the 
Manna model (Fig. 2), whose conditional distributions 
appear to be symmetric. Data from avalanche areas up 
to a ~ 10^ provide striking evidence that both models 
share the same asymptotic behavior with an exponent 
jsa = 1-35 ± 0.05, in agreement with other published re- 
sults P, |l0| , p^pT| . The scahng of the other expectation 
values is also biased as it is apparent from the bending 
in the curves reported in Refs. The correction 

of the bias is not so straightforward as in the case we 
have discussed but can be obtained from the analysis of 
p{x\y). This discussion clearly shows that conditional ex- 
pectation values are not a reliable method to determine 
the universality class of a model. 

To confirm the conclusion about a single universality 
class we perform the moment analysis on the distribu- 
tions P{x, L) in close analogy with the recent work of De 
Menech et al. on the two-dimensional BTW model. 
Here we apply the moments analysis on both BTW 
and Manna models, taking advantage of the large sizes 
reached in our numerical simulations. We define the q- 
moment of a; on a lattice of size L as {x'^)l — J x'^P(x)dx. 
If FSS hypothesis (Eq. (|^)) is valid, at least in the asymp- 
totic limit {x 00), we can transform z = x/L^'' and 
obtain 

= f z'^+''g{z)dz - (6) 



or in general {x'')l L'^^^*). The exponents ax{q) can 
be obtained as the slope of the log-log plot of < x'' >l 
versus L. Using Eq. (^), we obtain {x'^'^^) l / {x'^) l ^ L^'' 
or (Jx{q+ 1) — fyx{q) = Px, so that the slope of ax{q) as 
a function of q is the cutoff exponent Px = dax(q) / dq. 
This is in general not true for small q because the integral 
in Eq. (|^) is dominated by the lower cutoff. In particular, 
corrections to scaling of the type ~ L'^'-'^i^ F{L) are 

important for q < Tx — 1. For instance, when q ~ r^, — 1, 
logarithmic corrections give rise to effective exponents up 
to very large lattice sizes. Finally, normalization imposes 
<Jx{0) = 0. 

In Fig. 3 we show the results obtained from the mo- 
ment analysis of the distribution P{s) for the Manna and 
the BTW model. In this case we can use the exact re- 
sult (s) ~ L^, which implies (Ts(1) = 2, as a test for the 
convergence of our simulations to the asymptotic scaling 
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regime. This relation is fulfilled and the (Ts{q) of the two 
models are indistinguishable for g > 1, indicating uni- 
versal scaling behavior. We observe small deviations for 
small q which are due to the non-universal lower cutoff. 
By measuring slope of as{q), we obtain /3s ~ 2.7. This 
value is larger than the value reported in Ref. ||l5| (i.e. 
D ~ 2.5), where small lattice sizes have been used. We 
have repeated the same analysis for the P{T, L) and the 
P(a, L) and the measured cutoff exponents fit and Pa are 
reported in Table I. Also in this case the exponents for 
the two models share the same values within error bars, 
confirming the presence of a single universality class. 

As a final consistency test, we used the data col- 
lapse method in order the check the FSS hypothesis, 
which states that rescaling qx = xjL^'^ and Pq^ = 
P{x, L)L^'='^^ , the data for different L must collapse onto 
universal curves. If FSS is verified, we can compute the 
exponent Tx from the scaling relation [2 — Tx)I3x = 
that should be satisfied for enough large sizes. Using 
the values of Px reported in Table I and the values ob- 
tained for cra;(l) we find the exponents Tx to be inserted 
in the data collapse. For instance, using the exact re- 
sult crs(l) = 2 and the estimated (3s = 2.7, we obtain 
Ts = 1.27. The data collapse with these values is satis- 
factory for both models (see Fig. 4). 

In the same way, we obtain very good data collapse for 
the Manna model P{a) and P{t) distributions, yielding 
Tf = 1.5 and = 1.35. On the other hand, we find that 
the BTW data collapses for time and area distributions 
are not compatible with the FSS hypothesis. The linear 
behavior of the moments analysis, however, ensures that 
for large sizes the FSS form must be approached. This 
result can be explained if we assume that the scaling in 
the BTW model displays subdominant corrections of the 
form P{x) = (Cix-^i + C2X-^^ + ...)g{x/xc), where 
are non-universal constants. These corrections are com- 
patible with the linear behavior at large q, but the decay 
of the P{x) is not a simple power law for small x and 
thus FSS is not obeyed. It is worth to remark that the 
time and area distributions span over much less order 
of magnitude than the size distribution, which could ex- 
plain why subdominant corrections are more relevant in 
the first two cases. Subdominant corrections are due to 
higher order operators in the dynamics and do not deter- 
mine the universality class, since the asymptotic scaling 
behavior is ruled by the leading power. 

In summary, we have presented strong numerical evi- 
dence pointing towards a single universality class for the 
Manna and the BTW model. In particular, we show 
that previous analyses [|ll],0 are not reliable because of 
systematic biases introduced by the method employed. 
Further work is needed in order to quantify the extent of 
subdominant corrections to scaling in the BTW model. 
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Model 




Pt 


Pa 




Manna 
BTW 


2.74 ± 0.02 
2.73 ± 0.02 


1.50 ±0.02 
1.52 ±0.02 


2.02 ±0.02 
2.01 ±0.02 


1.27 ±0.01 
1.27 ±0.01 



TABLE I. Values of the critical exponents describing the 
scaling of the cutoff of the distributions for different models 
in d = 2. The results are obtained from the moments analysis 
(see text) . Note that the exponents /3s , /3t and (3a are usually 
reported in the literature as D, z and d/, respectively. 
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FIG. 1. The figure shows the probability distribution of 
having an avalanche size s given its area a for the BTW model. 
The inset shows that all data collapse onto the universal scal- 
ing function p{s\a) = a~^°"/((s — a))/a^°°), with jsa — 1.35. 



FIG. 3. Plot of as{q) for the BTW and Manna model. The 
linear part has slope 2.74. Note the non universal corrections 
to the linear behavior expected for g ~ r — 1 ~ 0.3. 
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FIG. 4. Data collapse analysis of the avalanche size distri- 
bution for the Manna and BTW (inset) models. The values 
used for the critical exponents are — 1.27 and /3s — 2.7. 
Lattice sizes used are reported in figure. 



FIG. 2. Conditional expectation value E{s\a) for the BTW 
and Manna model (after bias subtraction) . The slope is given 
by 7sa ~ 1-35 ± 0.05 for both curves. 
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